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SMALL GRAINS FEATURE SELECTION 


INTRODUCTION 


Separation of spring wheat (SW) from compe- 
ting small grains (SG) hrs been the subject of several 
studies. In this project we wish to approach the problem 
from an emperical point of view, using nonparametric 
discriminate function methods to investigate the feasi- 
bility of SW/SG separation. A limited- data set consis- 
ting of five segments with multiple acquisitions was 
used to illustrate the software and to make preliminary 
conclusions with respect to appropriate features. The 


data sets 

used are 

given 

in 

Table 3 




Segment 

1 

Acquisition data in 1976 
2 3 4 5 6 

7 

1614 

130 



183 

201 

219 


1618- 

127 


163 


199 


235 

1624 

128 

146 





236 

1642 

' 127 

145 

163 

182 

199 


236 

1645 

127 

145 

164 

181 



235 


Table 1. Acquisitions use in experiment. 
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LINEAR DISCRIMINATE FUNCTIONS 


Several methods including Principal Components , 

Fischer’s linear discriminate function, minimization of 

the Perceptron criteria function, and Minimum Squared 

Error (MSE) procedures were considered and tested. The 

MSE procedure using the Ho-Kashyap algorithm was chosen 

as the most consistent, based on limited testing. A brief 

discription of the procedure follows. 

* T 

Let ( *^*i 1 / 2^i2 / ... / be the 

measurement vector for the ith prototype. Define 


Y i 


( 1, xt ) 


( -1, -X* ) 


for X^ in class 1 
for X^ in class 2 


The problem is to select a (p+1) vector a , such that 



for all i ,= 1, 


n. 


Or, failing that, minimize the number of errors. A 
more tractable problem is the MSE criteria which is 
defined as follows. 



3 . 


Let 


A 




n x (p+1) 


Now find a solution (if possible) to the equation 

A a = 3 > 0 . 

Or, minimize '{ |Aa - 3| | over all a and 3 • For 
each fixed 3 , the minimum. norm solution is given by 

a = A* 3/ 

# 

where A is the pseudo- inverse of A . So the 
difficulty is one of determining the appropriate 3 . 
The Ho-Kashyap algorithm solves this problem in the 
following way. Define the initial values of two vector 
sequences a and 3 by 

3(0) = (1, 1) T 

= A # 3(0) . 


a(0) 
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For k = 1,2,.... define 

e (k) = A o(k) - $(k) 

e + {k) = J{ e (k) + je(k) | ) . 

and 

B (k+1) = B (k) + pe + (k) 

a (k+1) = a(k) + pA # e + (k) 

where p > 0. If the prototypes are separable, then 
| | e (k) | | converges to 0. If not, then | |e(k) | | converges 
to c > 0 . 

THE EXPERIMENT 


In this test we have used the Ho-Kashyap 
algorithm to determine three discriminant functions for 
each of five segments over various pass combinations. 

The three discriminant functions are defined as follows. 
4-CH The original four channel LANDSAT measurements 

are used for each acquisition. 

L (B, G) The Tassel Cap coordinates B and G are 
computed for each acquisition. 

Q (B,G) The Tassel Cap coordinates B and G plus the 

2 2 

quadratic terms B , G , and BG are computed for each 
acquisition. In Tables 2-6 the error rates are given for 
each of the above discriminant f unctions , calculated by 
the Ho-Kashyap algorithm. 
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SEGMENT 1614 
ERRORS (SW/SG) 


ACQ. 

4-CH 

L (B , G) 

Q(B,G) 



1 

1/5 

1/11 

1/11 

4 

3/7 

2/7 

1/7 

5 

2/5 

1/6 

0/6 

6 

2/3 

0/6 

0/4 

1,4 

1/4 

1/6 

1/6 

1,5 

0/0 

1/5 

0/5 

1,6 

0/1 

2/5 

0/4 

4 , 5 

3/4 

1/5 

2/5 

4,6 

2/3 

1/4 

0/2 

5,6 

1/3 

0/4 

0/1 

1,4,5 

0/0 

1/5 

1/2 

1,4,6 

0/0 

0/3 

0/2 

1,5,6 

0/0 

0/5 

0/2 

4,5,6 

2/3 

1/5 

o/i 

1,4,5, 6 

0/0 

2/4 

0/0 


Labeled Dot 

Distribution: 

SW-31 / 

SG-13 

Acquisition 

Dates: 

1 

76130 



4 

76183 



5 

76201 



6 

76219 


Table 2. 
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SEGMENT 1618 


ERRORS (SW/SG) 


ACQ. 

4-CH 

L(B,G) 

Q(B,G) 

1 

0/22 

0/22 

2/21 

3 

3/21 

2/21 

4/14 

5 

3/10 

3/15 

3/4 

7 

4/9 

2/13 

2/13 

1/3 

2/20 

2/20 

6/13 

1/5 

3/10 

3/15 

2/3 

1/7 

5/9 

2/12 

4/11 

3/5 

3/10 

3/15 

1/4- 

3/7 

4/7 

1/13 

2/7 

5/7 

3/8 

3/13 

1/5 

1/3/5 

3/9 

3/13 

1/3 

1/3/7 

5/7 

4/11 

5/5 

1,5/7 

4/5 

3/9 

0/2 

3/5,7 

3/4 

3/12 

1/2 

1/3,5, 7 

3/3 

3/9 

0/2 


Labeled Dot Distribution: SW-60 / SG-22 


Acquisition Dates: 1 - 76127 

3 - 76163 
5 - 76199 
7 - 76235 


Table 3 
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SEGMENT 1624 
ERRORS (SW/SG) 


ACQ. 

4-CH 

L(B,G) 

Q (B/G) 

1 

4/18 

0/22 

1/20 

2 

0/22 

0/21 

1/21 

7 

5/16 

5/16 

5/13 

1/2 

3/16 

0/21 

1/19 

1/7 

6/11 

5/14 

4/13 

2/7 

5/13 

5/13 

5/13 

1/2,7 

4/12 

4/14 

4/11 


Labeled Dot Distribution: SW-62 / SG-22 

Acquisition Dates: 1 - 76128 

2 - 76146 

7 - 76236 


Table 4. 
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SEGMENT 1642 
ERRORS (SW/SG) 


ACQ 

4-CH 

L(B,G) 

Q(b,g) 

1 

0/18 

0/19 

0/18 

2 

2/14 

2/18 

1/13 

3 

3/16 

2/15 

3/12 

4 

0/19 

0/19 

0/16 

5 

0/18 

0/19 

3/16 

7 

1/19 . 

1/18 

0/18 

1/2 

3/14 

1/18 

1/12 

1,3 

2/9 

2/12 

4/11 

1,4 

1/17 

0/19 

0/16 

1,5 

3/16 

2/18 

4/9 

1,7 

1/17 

1/17 

1/16 

2,3 

3/11 

3/12 

1/9 

2,4 

3/18 

2/18 

2/11 

2,5 

2/14 

2/18 

5/10 

2,7 

4/15 

1/18 

1/10 

3,4 

3/13 

3/14 

3/10 

3,5 

5/10 

3/14 

4/8 

3,7 

4/13 

4/14- 

2/8 

4,5 

2/16 

0/18 

2/19 

4,7 

1/19 

1/18 

0/15 

5,7 

0/17 

0/18 

3/14 

1,2, 3, 4, 5, 7 

7/7 

4/9 

0/0 

Labeled Dot 

Distribution : 

SW-58 / 

SG-19 

Acquisition 

Dates: 

1 

76127 



2 

76145 



3 

76163 



4 

76182 



5 

• 76199 



7 

76236 


Table 5 



9 . 


SEGMENT 1645 
ERRORS (SW/SG) 


ACQ 

4-CH 

L(B,G) 

Q (B,G) 

1 

0/20 

0/20 

0/20 

2 

0/20 

0/20 

0/20 

3 

1/20 

0/20 

0/19 

4 

1/19 

0/20 

0/20 

7 

0/19 

0/20 

2/16 

1,2 

0/19 

0/20 

0/20 

1,3 

0/18 

' 0/20 

1/19 

1,4 

0/17 

0/20 

0/20 

1,7 

1/12 

2/17 

2/11 

2,3 

0/20 

0/20 

0/19 

2,4 

1/18 

0/20 

0/20 

. 2,7, 

0/18 

0/20 

3/11 

3,4 

4/16 

2/20 

1/18 

3,7 

2/16 

1/17 

3/11 

4,7 

3/10 

3/13 

3/10 

1/ 2 , 3, 4 , 7 

3/10 

6/11 

2/6 

Labeled Dot 

Distribution: 

SW-75 / 

SG-20 

Acquisition 

Dates: 

1 ' - 

76127 



2 

76145 



3 

76164 



4 

76181 



7 

76235 


Table 6 
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CONCLUSIONS 


The first observation must be, that, in these 
segments, the separation of SW from SG is not an easy 
task. A reasonable error rate was achieved with one 
pass only in segment 1614. Reasonable two pass error 
rates were achieved in segments 1614 and 1618, and to a 
lesser degree in 1642. Segment 1624 had inadequate 
acquisitions and segment 1645 gave poor results even 
when all acquisitions were used. Segment 1645 and 1642 
did not have acquisitions in windows 5 or 6 . These two 
windows provided the best results in the other three 
segments in single pass or two pass combinations. 

The other observation is that generally the . 
Q(B,G) features provided as good or better separation 
as did the 4-CH features. The L (B, G) • features did not 
compete as well. The advantage of the Q(B,G) features 
is that the two dimensional quadratic discriminant 
function 'can be plotted on a per pass basis, making the 
prospect of generating graphical AI aids a possibility. 

In conclusion, it appears that windows 5 and 6 
play an important role in the SW/SG separation problem. 
(The corresponding crop indices should be determined 
for this strata.) In addition we recommend further 
testing of the Q(B,G) features over a broad range of 
spring wheat blind sites. 
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**** PROGRAM UHLDF **** 

THIS PROGRAM COMPUTES THE 2-CLASS NCNPARAMETRIC 
DISCRIMINATE FUNCTION. THE HO-KASHYAP ALGORITHM 
IS IMPLEMENTED. 


DIMENSION D (4000) ,DINV(4000) ,A (40 ), ICLASS (500) ,B(500) ,Y(500) 
DIMENSION Z (500) ,C (500) ,IQ (500) ,V (40) ,E (500) 

DIMENSION U (40,40) ,AFIAG(40) f ATEMP (40) ,W(25D0) 

EQUIVALENCE (U (1,1) ,B (1) ) , (AFIAG (1) ,Z (1) ) 

* , (ATEMP (1) ,Y (51) ) , (W(1),B(1)) 

DATA YES/’Y'/ 

66 WRITE (108 ,6600) 

READ (105 ,6666) YESS 
5600 FORMAT {’ AGAIN?’) 

5666 FORMAT (Al) 

IF (YESS. NE. YES) STOP 

GET DATA ARRAY 

D - NSAMP BY NV ARRAY CONTAINING PROTOTYPES AS RONS. 
DINV - TRANSPOSE OF THE PSEUDO-INVERSE OF D. (NSAMP BY 

REWIND 1 

CALL GETD (DINV, NSAMP ,NV, ICLASS) 

IF(NSAMP.EQ.O) GO TO 500 
NSIZE=NSAMP*NV 
OUTPUT, NSIZE 

CALL TRANS (DINV,D,NV,NSAMP) 

CALL MOVE (D, DINV, NSIZE) 


MAXITR=100 


COMPUTE PSEUDO-INVERSE OF D. (TRANSPOSE) 

INITIALIZE B— VECTOR 

INITIALIZE A-VECTOR. (LINEAR DISCRIMINATE FUNCTION) 
COMPUTE MISSCLASSIFICATIONS FOR INITIAL DF. . 

CALL GINV2M (DINV ,NSAMP , NV,NSAMP , NV , KZ , U , AFIAG , ATEMP , 1 . E-l 2 ) 
CALL FILL (B, NSAMP, 1.) 

CALL TPRD (DINV, B,A,NSAMP,NV, 0,0,1) 

WRITE (108, 8300) (A(J) ,J=1,NV) 

CALL MPRD (D, A,Y, NSAMP ,NV, 0, 0 ,1) 

CALL MISSCL (Y ,NSAMP ,MISSl ,MISS2 , ICLASS) 
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MISS=MISS1+MISS2 

UHLD0500 


OUTPUT MISS1,MISS2,MISS 

UHLD0510 

c 


UHLD0520 

c 

CALL HO-KASHYAP ALGORITHM. 

UHLD0530 

c 


UHLD0540 


CALL HOKASH (D , DINV r A, B ,Y , E ,NV,NSAMP ,MAXITR, ICLASS) 

UHLD0550 

8800 FORMAT (2X,5FlO. 4) 

UHLD0560 


GO TO 66 

UHLD0570 


END 

UHLD0580 

C 


UHLD0590 

C 

COPY ARRAY X INTO ARRAY Y 

UHLD0600 

C 


UHLD0610 

C 


UHLD0620 


SUBROUTINE MOVE (X,Y,N) 

UHLD0630 


DIMENSION X{1) ,Y(1) 

UHLD0640 


DO 1 1=1, N 

UHLD0650 


1 Y(I)=X{I) 

DHLD0660 


RETURN 

UHLD0670 


END 

UHLD0680 

c 


UHLDO'690 

c 

INITIALIZE N LOCATIONS IN ARRAY X WITH VALUE C 

UHLD0700 

c 


UHLD0710 

c 


UHLD0720 


SUBROUTINE FILL(X,N,C) 

UHLD0730 


DIMENSION X(l) 

UHLD0740 


DO 1 1=1, N 

UHLD0750 


1 X(I)=C 

UHLD0760 


RETURN 

UHLD0770 


END 

UHLD0780 

c 


UHLD0790 

c 

MOVE TRANSPOSE OF MATRIX A INTO B. 

UHLD0800 

c 


UHLD0810 

c 


UHLD0820 

c 


UHLD0830 


SUBROUTINE TRANS (A ,B , N ,M) 

UHLD0840 


DIMENSION A(N,M) ,B (M,N) 

UHLD0850 


DO I 1=1, N 

UHLD0860 


DO I J=1,M 

UHLD0870 


1 B (J,I)=A(I ,J) 

UHLD0880 


RETURN 

UHLD0890 


END 

UHLD0900 
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HO-KASHYAP ALGORITHM. 


SUBROUTINE HOKASH (D ,DINV,A,B,Y,E ,NV,NSAMP ,MAXITR, ICLASS) 
DIMENSION D (1) ,DINV(1) ,A(1) ,B(1) ,Y(1) ,E(1) , ICLASS (1) 
ITR-0 

10 ITR=ITR+1 

CALL MPRD(D,A,Y,NSAMP,NV, 0,0,1) 

CALL MSUB(Y,B,E,NSAMP, 1,0,0) 

CALL MISSCL (Y ,NSAMP ,MISSl f MISS2 , ICLASS ) 

IF (ITR.GT .MAXITR) GO TO 200 
CALL TEST (E,NSAMP, KEY) 

IF (KEY) ^00,100,100 


ONE MORE TIME 


ORIGINAL PAGE IS 
OF POOR QUALITY 


100 CALL POS (E, NS AMP) 

CALL MAH) (B,E,B,NSAMP, 1,0,0) 

CALL TPRD (DINV,E,Y,NSAMP,NV, 0,0,1) 
CALL MADD(A,Y,A,SV,1,0,0) 

GO TO 10 


TERMINATE 


200 CONTINUE 

OUTPUT, ITR,MISS1,MISS2 
WRITE (108 , 8800 ) (A (J) , J=1 ,NV) 
1800 FORMAT {2X,5F10. 4) 

RETURN 

END 


SUBROUTINE TEST (Y,N, KEY) 
DIMENSION Y (N) 

KEY-0 

LN=0 

LP=0 

CALL NORM (Y ,N ,YNORM) 

DO 10 1=1, N 

IF (Y (I)/YNORM-1.E-50) 1,1,2 

1 IN=1 
GO TO 3 

2 LP=1 

3 IF (LN.EQ.l.AND.LP.EQ.l) RETURN 


HOKA0000 

HOKA0010 

HOKA0020 

HOKA0030 

HOKA0040 

HOKA0050 

HOKA0060 

HOKA0070 

HOKA0080 

HOKA0090 

HOKA0100 

HOKA0110 

HOKA0120 

HOKA0130 

HOKA0140 

HOK40150 

HOKA0160 

HOKA0170 

HOKA0180 

HOKA0190 

HOKA0200 

HOKA0210 

HOKA0220 

HOKA0230 

HOKA0240 

HOKA0250 

HOKA0260 

HOKA0270 

HOKA0280 

HOKA0290 

HOKA0300 

HOKA0310 

HOKA0320 

HOKA0330 

HOKA0340 

HOKA0350 

HOKA0360 

HOKA0370 

HOKA0380 

HOKA0390 

HOKA0400 

HOKA0410 

HOKA0420 

HOKA0430 

HOKA0440 

HOKA0450 

HOKA0460 

HOKA0470 

HOKA0480 

HOKA0490 
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10 CONTINUE 

HOKA0500 


KEY=1 

HOKA0510 


IP (LP.EQ.O) KEY=-1 

HOKA0520 


RETURN 

HOKA0530 


END 

HOKA0540 

c 


HOKA0550 

c 


HOKA0560 

c 


HOKA0570 

c 


HOKA05S0 

c 


HOKA0590 


SUBROUTINE POS(Y,N) 

HOKA0500 


DIMENSION Y (N) 

HOKA0610 


DO 10 1=1 ,N 

HOKA0620 


10 Y{I)=(Y(I)+ABS(Y(I)))/2. 

HOKA0630 


RETURN 

HOKA0640 


END 

HOKA0650 

c 


HOKA0660 

c 


HOKA0670 

c 


BOKA0680 

c 


HOKA0690 

c 


HOKA0700 


SUBROUTINE MISSCL (Y , NS AMP f MISSl ,MISS2 , ICIASS) 

HOKA0710 


DIMENSION Y(l) , ICIASS (1) 

HOKA0720 


MISS1=0 

HOKA0730 


MISS 2=0 

HOKAO740 


DO 30 1=1, NS AMP 

HOKA0750 


IF(Y(I).G5.0.) GO TO 30 

HOKA0760 


IF (ICIASS (I) .EQ.l) MISS 1 '"MISS 1+1 

HOKA0770 


IF (ICIASS (I ) .EQ. 2) MISS2=MISS2+1 

HOKA07S0 


30 CONTINUE 

HOKA0730 


RETURN 

HOKA0800 


END 

BOKA0810 

c 


HOKA0820 

c 


HOKA0830 

c 


HOKA0840 

c 


HOKA0850 

c 


HOKA0860 


SUBROUTINE NORM(Y ,N,YNORM) 

HOKA08“0 


DIMENSION Y(l) 

HOKA0880 


YNORM=0 

HOKA0890 


DO 10 1=1, N 

HOKA0900 


10 YNORM=YNORM+Y(I)*Y(I) 

HOKA0910 


YNORM=SQRT (YNORM) 

HOKA0920 


RETURN 

HOKA0930 


END 

HOKA0940 
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*** 4 CHANNEL LANDSAT *** 

THIS VERSION OF GETD USES THE 4 LANDSAT CHANNELS AS 
FEATURES. (UP TO 7 PASSES) 

D - THIS IS THE DATA ARRAY WHICH IS RETURNED TO THE 
CALLING ROUTINE. THE DIMENSION IS NV BY NSAMP. 

THE COLUMN VECTORS OF D CONSIST OF 

( 1. , XI, X2, ... , XN) 

IF THE PROTOTYPE IS IN CLASS 1 AND THE NEGATIVE 
OF THE ABOVE VECTOR IF THE PROTOTYPE IS IN 
CLASS 2. 

NSAMP- THE -NUMBER OF PROTOTYPES. THIS VALUE IS RETURNED TO 
THE CALLING PROGRAM. 

NV - THE NUMBER OF VARIABLES PLUS ONE FOR THE 

CONSTANT TERM. THIS VALUE IS RETURNED TO THE CALLING 
PROGRAM. 


GD010000 

GD010010 

GD0I0020 

GD010030 

GD010040 

GD010050 

GD010060 

GD010070 

GD010080 

GD010090 

GD010100 

GD010110 

GD010120 

GD010130 

GD010140 

GD010150 

GD010160 

GD010170 

GD010180 

GD010190 

GD010200 

GD010210 

GD010220 

GD010230 

GD010240 


ICLASS- THIS IS A VECTOR, RETURNED TO THE CALLING PROGRAM WHICH 
IDENTIFIES (1 OR 2) THE CLASS .ASSIGNMENT OF EACH 
PROTOTYPE. 


SUBROUTINE GETD (D , NSAMP ,NV , ICLASS ) 
DIMENSION D (1) , ICLASS (1) ,X (28) ,IP (7) 

DATA IBLK,IO,IW/' ' ,'0','V!'/ 

88 WRITE (108,4000) 

4000 FORMAT (' INPUT NUMBER OF PASSES. 1-7*) 
READ (105,3000 ,ERR=8 8) NPASS 
OUTPUT NPASS 
WRITE (108 ,2000) 

2000 FORMAT (' INPUT PASS NO. 1-7') 

READ (105,3000, ERR=88) (IP(KK) ,KK=1, NPASS) 
WRITE (108,2500) (IP (KK) ,KK=1, NPASS) 

2500 FORMAT (7X, 713) 

3000 FORMAT (711) 

NSAMP=0 

N1=0 

N2=0 

IND=0 

NV=NPASS*4+1 
DO 100 K=l,209 


GD010250 

GD010260 

GD010270 

GD010280 

GD010290 

03010300 

GD0103I0 

GD010320 

GD010330 

GD010340 

GD010350 

GDG10360 

GD010370 

GD010380 

GD010390 

GD010400 

GD0104IO 

GD010420 

GD010430 

GD010440 

GD010450 

GD010460 

GD010470 

GD010480 

GD010490 



non 
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READ LABELED DOTS (UNDER FORMAT 1000) 

READ (1, 1000, END=999) IG1,IG2,X 
1000 FORMAT (2Al,8F5.1,/,2X,8F5.1,/,2X,8F5.1,/ r 2X,4F5.1) 
IF(IGl.NE.IBLK) GO TO 100 
IF (IG2.EQ.IBLK.OR.IG2.EQ.IO) GO TO 100 
DO 10 KK=1,NPASS 
11= (IP(KK)-1)*4+1 
IF (X (11) .GT.98.5) GO TO 100 
10 CONTINUE 

NSAMP=NSAMP+1 

ICLASS (NSAMP )=1 

IF (IG2.NE.IW) ICLASS (NSAMP) =2 

IND=IND+1 

D(IND)=1. 

IF (ICLASS (NSAMP) .EQ. 2) D (IND)=-1 . 

DO 30 NPP=1 r NPASS 

11= (IP (NPP) -1 ) *4+1 

12=11+3 

DO 20 J=Il,l2 

IND=IND+1 

D(IND)=X(J) 

IF (ICLASS (NSAMP) .EQ. 2) D (IND)=-D (HO) 

20 CONTINUE 
30 CONTINUE 

IF (ICLASS (NSAMP) .EQ.l) NI=Nl+l 
IF (ICLASS (NSAMP) .EQ. 2) N2=N2+1 
100 CONTINUE 
999 CONTINUE 

OUTPUT N1.N2, NSAMP 

RETURN 

END 


GD010500 

GD010510 

GD010520 

GD010530 

GD010540 

GD010550 

GD010560 

GD010570 

GD010580 

GD010590 

GD010600 

GD010610 

GD010620 

GD010630 

GD010640 

GD010650 

GD010660 

GD010670 

GD010680 

GD010690 

GD010700 

GD010710 

GD010720 

GD010730 

GD010740 

GD010750 

GD010760 

GD010770 

GD010780 

GD010790 

GD010800 

GD010810 

GD010820 
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C GD020000 

C GD020010 

C *** B,G *** GD020020 

C - GD02G030 

C THIS VERSION OF GETD USES THE FIRST 2 TASSEL CAP VARIABLES (B,G) GD020040 

C FEATURES. (UP TO 7 PASSES) GD020050 

C GD020060 

C D - THIS IS THE DATA ARRAY WHICH IS RETURNED TO THE GD020070 

C CALLING ROUTINE. THE DIMENSION IS NV BY NSAMP. GD020080 

C GD020090 

C THE COLUMN VECTORS OF D CONSIST OF • GD020100 

C GD020110 

C ( 1. , XI, X2, ... , XN) GD020120 

C GD020130 

C IF THE PROTOTYPE IS IN CLASS 1 AND THE NEGATIVE GD020140 

C OF THE ABOVE VECTOR IF THE PROTOTYPE IS IN GD020150 

C CLASS 2. GD020160 

C GD020170 

C NSAMP- THE NUMBER OF PROTOTYPES. THIS VALUE IS RETURNED TO GD020180 

C THE CALLING PROGRAM. GD020190 

C GD020200 

C NV - THE NUMBER OF VARIABLES PLUS CHE FOR THE GD020210 

. C CONSTANT TERM. THIS VALUE IS RETURNED TO THE CALLING GD020220 

C PROGRAM. GD020230 

C • GD020240 

C ICLASS- THIS IS A VECTOR, RETURNED TO THE CALLING PROGRAM WHICH GDQ20250 

C IDENTIFIES (1 OR 2) THE CLASS ASSIGNMENT OF EACH GD020260 

C PROTOTYPE. GD020270 

C GD020280 

C ' GD020290 

C GD020300 

SUBROUTINE GETD (D ,NSAMP ,NV, ICLASS) GD020310 

DIMENSION D (1) , ICLASS (1) ,X(28) ,IP(7) GD020320 

DATA IBLK,IO,IW/* VOW/ GD020330 

88 WRITE (108,4000.) GD020340 

4000 FORMAT (* INPUT NUMBER OF PASSES. 1-7') GD020350 

READ (105,3000 ,ERR=88) NPASS GD020360 

OUTPUT NPASS GD020370 

WRITE (108,2000) GD020380 

2000 FORMAT (' INPUT PASS NO. 1-0 ',) GD020390 

READ (105 ,3000, ERR=88) (IP(KK) ,KK=1, NPASS) GD020400 

WRITE (108,2500) (IP (KK) ,KK=1 , NPASS-) GD020410 

2500 FORMAT (7X, 713) GD020420 

3000 FORMAT (711) GD020430 

NSAMP=0 GD020440 

N1=0 GD020450 

N2=0 GD020460 

IND=0 GD0204O0 

NV=NPASS*2+1 GD020480 

DO 100 K=l,209 GD020500 

C GD020510 
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READ LABELED DOTS (UNDER FORMAT 1000) 

READ (1, 1000, END=999) IG1,IG2,X 
1000 FORMAT (2Al,8F5.1,/,2X,8F5.1,/,2X,8F5.1,/,2X,4F5.1) 
IF (IG1.NE.IBLK) GO TO 100 
IF (IG2.EQ.IBLK.OR.IG2.EQ. IO) GOTO 100 
DO 10 KK=1 ,NPA3S 
Il= (IP (KK) -1 ) *4+1 
IF (X (II) .GT. 98. 5) GO TO 100 
10 CONTINUE 

NSAMP=NSAMP+1 

ICLASS (NS AMP) -1 

IF (IG2.NE.IW) ICIASS (NSAMP) =2 

IND=IND+1 

D(IND)=1. 

IF (ICIASS (NSAMP) .EQ.2) D(IND)=-1. 

DO 30 NPP=1,NPASS. 

11= (IP (NPP) -1 ) *4+1 
CALL KAUTH (X (11) ,B,G) 

D(IND+1)=B 
D (IND+2)=G 
DO 20 J=l,2 

IF (ICIASS (NSAMP) .EQ.2) D (IND+J)=~D (IRD+J) 

20 CONTINUE 
IND=IND+2 
30 CONTINUE 

IF (ICLASS (NSAMP) .EQ.l) Nl=$!l+1 
IF (ICIASS (NSAMP) .EQ.2) N2=N2+1 
100 CONTINUE 
999 CONTINUE 

OUTPUT Nl,N2 f NSAMP 

RETURN 

END 


GD020520 

(33020530 

GD020540 

GD020550 

GD020560 

GD020570 

GD020580 

GD020590 

GD020600 

GD020S10 

GD020620 

GD020630 

GD020640 

GD020650 

GD020660 

GD0206'70 

GD020680 

GD020690 

GD0 20700 

GD020710 

GD020720- 

GD020730 

GD020740 

GD020750 

GD020760 

GD0207*70 

GD020780 

GD020790 

GD020800 

GD020810 

GD020820 

GD020830 

GD020840 
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**** B,G QUAD *** 

THIS VERSION OF GETD USES THE FIRST 2 TASSEL CAP VARIABLES PLUS 
QUADRATIC TERMS AS FEATURES.' (UP TO 7 PASSES) 

D - THIS IS THE DATA ARRAY WHICH IS RETURNED TO THE 
CALLING ROUTINE. THE DIMENSION IS NV BY NSAMP . 

THE COLUMN VECTORS OF D CONSIST OF 

( 1. , XI, X2, ... , XN) 

IF THE PROTOTYPE IS IN CLASS 1 AND THE NEGATIVE 
OF THE ABOVE VECTOR IF THE PROTOTYPE IS IN 
CLASS 2. 

NSAMP- THE NUMBER OF PROTOTYPES. THIS VALUE IS RETURNED TO 
THE CALLING PROGRAM. 

NV - THE NUMBER OF VARIABLES PLUS ONE FOR THE 

CONSTANT TERM. THIS VALUE IS RETURNED TO THE CALLING 
PROGRAM. 

ICLASS- THIS IS A VECTOR, RETURNED TO THE CALLING PROGRAM WHICH 
IDENTIFIES (1 OR 2) THE CLASS ASSIGNMENT OF EACH 
PROTOTYPE. 


SUBROUTINE GETD (D ,NSAMP ,NV , ICLASS ) 
DIMENSION D(l) , ICLASS (1) ,X{28) ,IP(7) 

DATA IBLK,IO,IW/' 1 , 'O' , 'W'/ 

WRITE (108, 4000) 

FORMAT (' INPUT NUMBER OF PASSES. 1-7’) - 

READ (105, 3000, ERR=88) NPASS 
OUTPUT' NPASS 
WRITE (108,2000) 

FORMAT (' INPUT PASS NO. 1-7’) 
READ(105,3000,ERR=88) (IP(KK) ,KK=1, NPASS) 
WRITE (108,2500) (IP (KK) ,KK=1, NPASS) 

FORMAT (7X, 71 3) 

FORMAT (711) 

NSAMP=0 

N1=0 

N2=0 

IND-0 

NV=NPASS *5+1 
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DO 100 K=l f 209 


READ LABELED DOTS (UNDER FORMAT 1000) 

READ (1, 1000, END=999) IG1,IG2,X 
1000 FORMAT (2A1,8F5.1,/,2X,8F5.1,/,2X,8F5.1,/,2X,4F5.1) 
IF (IGl.NE.XBUC) GO TO 100 
IF (IG2.EQ.IBIK.OR.IG2.EQ.IO) GO TO 100 
DO 10 KK=1,NPASS 


11= (IP(KK)-1)*4+1 
IF (X (11) .GT.98.5) GO TO 100 
10 CONTINUE 

NSAMP=NSAMP+1 

ICLASS (NSAMP) =1 

IF(IG2.NE.IW) ICLASS (NSAMP ) ~2 

IND=IND+1 

D(IND)=1. 

IF (ICIASS (NSAMP) .EQ.2) D(IND)=-1. 
DO 30 NPP=1,NPASS 
11= (IP (NPP) -1) *4+1 
CALL KAUTH (X (II ) ,B ,G) 

D (IND+1)=B 
D (IND+2)=G 
D (IND+3 ) =B*B 
D(IND+4)=G*G 
D (IND+5 ) =B*G 

IF (ICLASS (NSAMP) .EQ.l) GO TO 25 
DO 20 J=IND+1, IND+5 
20 D (J)=-D (J) 

25 IND=IND+5 
30 CONTINUE 

IF (ICLASS (NSAMP) .EQ.l) Nl=Nl+l 
IF (ICIASS (NSAMP) .EQ. 2) N2=N2+1 
30 CONTINUE 


: page IS 


GD035000 
GD030510 
GD030520 
GD030530 
GD030540 
GD030550 
GD030560 
GD030570 
GD030580 
GD030590 
GD030600 
GD030610 
GD030620 
GD030630 
GD030640 
GD030650 
GD030660 
GD030670 
GDO 30680 
GD030690 
GD030700 
GD030710 
GD030720 
GD030730 
GD030740 
GD030750 
GD030760 
GD030770 
GD030780 
GD030790 
GD030800 
GD030810 
GD030820 
GD030830 


999 CONTINUE 

OUTPUT N1,N2,NSAMP 

RETURN 

END 


GDO 30 840 
GD030850 
GD030860 
GD030870 
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KAUTOOOO 

KAUTOOIO 

THIS ROUTINE COMPUTES THE TASSEL CAP COORDINATES KAOT0020 

3 AND G FROM THE FOUR LANDSAT CHANNELS IN THE VECTOR X. KAOT0030 

KAOT0040 

KAOT0050 

KAUT0060 

SUBROUTINE KAOTH (X,B,G) KAUT0070 

DIMENSION X (1) ,FB (4 ) ,FG {4) KAUT0080 

DATA FB, PG/. 33231,. 60316,. 67581 ,.26278, KAUT0090 

* -.28317 ,-.66006, .57735, ,38833/ KAUT0100 

3=0. KAUT0110 

G=0. KAOT0I20 

DO 1 1=1,4 KAOT0I30 

3=3+X ( I ) *FB ( I ) KAUT01 40 

1 G=G+X(I)*FG(I) KAUT0150 

RETURN KAUT0160 

END KAOT0170 
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SUBROUTINE GINV2M (A, MR , MC ,NR , NC ,KZ ,U , AFLAG, ATEMP ,TOK) GINVOOOO 

C— GINV0010 

C GINV0020 

C THIS ROUTINE COMPUTES THE TRANSPOSE OF THE GENERALIZED INVERSE GINV0030 
C OF A AND STORES IT IN A GINV0040 

C MR IS THE MAXIMUM ROW DIMENSION OF A GINV0050 

C MC IS THE MAXIMUM COLUMN DIMENSION OF A GINV0060 

C NR IS THE NUMBER OF ROWS IN A GINV0070 

C NC IS THE NUMBER OF COLUMNS IN A GINV0O8O 

C KZ OUTPUT, KZ= 0 IMPLIES THAT A IS SINGULAR GINV0090 

C KZ = 1 IMPLIES TEAT A IS NONSINGULAR GINV0100 

C U TEMPORARY WORKING STORAGE, DIMENSIONED MC BY MC GINV0110 

C AFLAG TEMPORARY WORKING STORAGE, DIMENSIONED MC GINV0120 

C ATEMP TEMPORARY WORKING STORAGE, DIMENSIONED MC GINV0130 

C TOK TOLERANCE FOR INNER PRODUCT ZERO GINV0I40 

C TOL TOLREANCE FOR THE RATIO OF TWO INNER PRODUCTS GINV0150 

C GINV0160 

C GINV0170 

DOUBLE PRECISION SUM , FAC GINV0180 

DIMENSION A (MR,MC) ,U (MC,MC) , AFLAG (MC) , ATEMP (MC) GINV0190 

KZ = 1 GINV0200 

DO 10 1=1, NO GINV0210 

DO 5 J=1,-NC GINV0220 

5 U (I ,J) =0.0 GINV0230 

10 0 (1,1) = 1.0 GINV0240 

DO 12 L=1,NC GINV0250 

FAC = DOT (MR,NR, A,L,L) GINV0260 

K=L GINV0270 

IF (FAC. GT. TDK) GO TO 11 GINV0280 

DO 13 J=1,NR GINV0290 

13 A (J ,L)=0.0 GINV0300 

KZ=0 GINV0310 

12 AFLAG (L)=0.0 GINV0320 

GO TO 1000 GINV0330 

11 FAC= 1.0/SQRT (FAC) GINV0340 

L=K GINV0350 

AFLAG (L ) =1 . 0 GINV0 360 

DO 15 1=1, NR GINV0370 

15 A(I,L) = A(I,L)*FAC GINV0380 

DO 20 1=1, NC GTNV0390 

20 U (I ,L) ' = U (I ,L) *FAC GINV0400 

C GINV0410 

C DEPENDENT COL TOLERANCE FOR N BIT FLOATING POINT FRACTION GINV0420 

C GINV0430 

N=2I GINV0440 

TOL=(10.*0.5**N)**2 GINV0450 

Ll=L+l GINV0460 

IF (Ll.LE.NC) GO TO 21 GINV047O 

DO 22 1=1, NR GINV04 80 

22 A (I ,L) = A (I ,L) *FAC GINV0490 
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GO TO 1000 


GINV0500 

21 DO 100 J=2,N0 

ORIGINAL: PAGE IS 

GINV0510 

DOTl = DOT {MR, NR/A, J,J) 
JMl = J-l 

£f©OI? QUALITY 

GINV0520 

GINV0530 

DO 50 L=l,2 


GINV0540 

DO 30 K=1,JM1 


GINV0550 

30 ATEMP{K) = DOT (MR,NR,A,J,K) 


GINV0560 

DO 45 K=l,JMl 


GINV0570 

DO 35 1=1, NR 


GINV0580 

35 A (I ,J) = A (I ,J) -ATEMP (K) *A(I ,K) *AFIAG (K) 

GINV0590 

DO 40 1=1, NO 


GINV0600 

40 U(I,J) = U(I,J)-ATEMP(K)*U(I,K) 


GINV0610 

45 CONTINUE 


GINV0620 

50 CONTINUE 


GINV0630 

DOT2 = DOT (MR,NR,A,J ,J) 


GINV0640 

IF ( (DOT2/DCT1) -TOL) 55,55,70 


GINV0650 

55 DO 61 1=1, JMl 


GINV0660 

SUM=0.0 


GINV06 7 0 

DO 60 K=l, I 


GINV0680 

60 SUM = SUM + U(K,I)*U{K,J) 


GINV0690 

61 ATEMP (I) = SUM 


GINV0 7 00 

DO 63 1=1, NR 


GINV0710 

SUM = 0.0 


GINV0 7 20 

DO 65 K=1,JM1 


GINV0730 

65 SUM = SUM - A(I,K) * ATEMP (K ) * AFLAC (K ) 

GINV0740 

63 A (I ,J) = SUM 


GINV0 7 50 

AFLAG(J) = 0.0 


GINV0 7 60 

KZ = 0 


GINV0770 

FAC = DOT (MC ,NC ,U , J , J) 


GINV0 7 80 

IF(FAC.GT.TOK) GO TO 66 


GINVO^O 

DO 67 1=1, NO 


GINV0800 

67 U(I,J) = 0.0 


GINV0810 

GO TO 100 


GINV0820 

66 FAC = 1.0/SQRT (FAC) 


GINV0830 

GO TO 75 


GINV0840 

70 IF (DOT2.GT.TOK) GO TO 71 


GINV0850 

KZ = 0 


GINV0860 

AFLAG(J) = 0.0 


GINV0870 

DO 72 1=1, NR 


GINV0880 

72 A (I ,J) = 0.0 


GINV0390 

GO TO 100 


GINV0900 

71 AFLAG(J) = 1.0 


GINV0910 

FAC = 1.0/SQRT (DOT2) 


GINV0920 

75 DO 80 1=1, NR 


GINV0930 

80 A (I ,J)=A (I , J) *FAC 


GINV0940 

DO 85 1=1, NO 


GINV0950 

85 U(I,J) = U { I , J ) *FAC 


GINV0960 

100 CONTINUE 


GINV09 7 0 

DO 130 J=1,NC 


GINV0980 

DO 130 1=1, NR 


GINV0990 
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FAC = 0.0 
DO 120 K=J,NC 

120 FAC = FAC+A (I ,K) *U (J ,K) 

130 A (I r J) = FAC 
1000 RETURN 
END 

FUNCTION DOT(MR,NR,A,JC,KC) 

C 

C COMPUTES THE INNER PRODUCT OF COLUMNS JC AND KC 
C OF MATRIX A 

DIMENSION A(MR f l) 

DOT = 0.0 
DO 5 1=1, NR 

5 DOT = DOT+ A (I , JC) *A (I ,KC) 

RETURN 

END 


GINV1000 

ginvioio 

GINV1020 

GINV1030 

GINVI040 

GINV1050 

GINV1060 

GINV1070 

GINV1080 

GINV1090 

GINV1100 

GINV1110 

GINV1120 

GTNV1130 

GINV1140 

GINV1150 
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c 



MADD 

10 

c 



• .MADD 

20 

c 



MACD 

30 

c 

SUBROUTINE MADD 


MADD 

40 

c 



MADD 

50 

c 

PURPOSE 


MADD 

60 

c 

ADD TWO MATRICES ELEMENT BY ELEMENT TO FORM 

RESULTANT 

MADD 

70 

c 

MATRIX 


MADD 

80 

c 



MADD 

90 

c 

USAGE 


MADD 

100 

c 

CALL MADD (A,B,R,N,M, MSA, MSB) 


MADD 

no 

c 



MADD 

120 

c 

DESCRIPTION OF PARAMETERS 


MADD 

130 

c 

A - NAME OF INPUT MATRIX 


MADD 

140 

c 

B - NAME OF INPUT MATRIX 


MADD 

150 

c 

R - NAME OF OUTPUT MATRIX 


MADD 

160 

c 

N - NUMBER OF ROWS IN A,B,R 


MADD 

170 

c 

M - NUMBER OF COLUMNS IN A,B,R 


MADD 

180 

c 

MSA - ONE DIGIT NUMBER FOR STORAGE MODE OF 

MATRIX A 

MADD 

190 

c 

0 - GENERAL 


MADD 

200 

c 

1 - SYMMETRIC 


MADD 

210 

c 

2 - DIAGONAL 


MADD 

220 

c 

MSB - SAME AS MSA EXCEPT FOR MATRIX B 


MADD 

230 

c 



MADD 

240 

c 

REMARKS 


MADD 

250 

c 

NONE 


MADD 

260 

c 



MADD 

270 

c 

SUBROUTINES AND FUNCTION SUBPROGRAM REQUIRED 


MADD 

280 

c 

IOC 


MADD 

290 

c 



MADD 

300 

c 

METHOD 


MADD 

310 

c 

STORAGE MODE OF OUTPUT MATRIX IS FIRST DETERMINED. ADDITION MADD 

320 

c 

OF CORRESPONDING ELEMENTS IS THEN PERFORMED. 

MADD 

330 

c 

THE FOLLOWING TABLE SHOTS THE STORAGE MOTE 

OF THE OUTPUT 

MADD 

340 

c 

MATRIX FOR ALL COl©INATIONS OF INPUT MATRICES 

MADD 

350 

c 

A 3 

K 

MADD 

360 

c 

GENERAL ' GENERAL 

GENERAL 

MATO 

370 

c 

GENERAL SYMMETRIC 

GENERAL 

MADD 

380 

c 

GENERAL DIAGONAL 

GENERAL 

MADD 

390 

c 

SYMMETRIC GENERAL 

GENERAL 

MADD 

400 

c 

SYMMETRIC SYMMETRIC 

SYMMETRIC 

MADD 

410 

c 

SYMMETRIC DIAGONAL 

SYMMETRIC 

MADD 

420 

c 

DIAGONAL GENERAL 

GENERAL 

MADD 430 

r 

DIAGONAL SYMMETRIC 

SYMMETRIC 

MADD 

440 

c 

DIAGONAL DIAGONAL 

DIAGONAL 

MADD 

450 

c 



MADD 

460 

c 

• 


. . .MADD 

470 

c 



MADD 

480 


SUBROUTINE MADD (A,B,R,N,M,MSA,MSB) 


MADD 

490 


DIMENSION A (1) ,B(1) ,R(I) 


MADD 

500 



non nan ono 


ORIGINAL PAGE IS 
Op POOR QUALITY 



MADD 510 

DETERMINE STORAGE MODE OF OUTPUT MATRIX 

MADD 520 


MADD 530 

IF (MSA-MSB) 1 , 5,1 

MADD 540 

5 CALL LDC(N,M,NM,N,M,MSA) 

MADD 550 

GO TO 100 

MAID 560 

7 MTEST=MSA*MSB 

MADD 570 

MSR=0 

MADD 580 

IF(MTEST) 20,20,10 

MADD 590 

10 MSR=T 

MADD 600 

20 IF (MTEST-2) 35,35,30 

MADD 610 

30 MSR=2 

MADD 620 


MADD 630 

LOCATE ELEMENTS AND PERFORM ADDITION 

MADD 640 


MADD 650 

35 DO 90 J— 1,M 

MADD 660 

DO 90 1=1, N 

MADD 670 

CALL LOC (I ,J,IJR,N,M,MSR) 

MADD 680 

IF (IJR) 40,90,40 

MADD 690 

40 CALL LOC (I ,J,IJA,N,M,MSA) 

MADD 700 

AEL=0 . 0 

MADD 710 

IF(IJA) 50,60,50 

MADD 720 

50 AEL=A (IJA) 

MADD 730 

60 CALL LOC (I ,J,IJ3,N,M,MSB) 

MADD "740 

BEL=0 . 0 

MADD 750 

IF(IJB) 70,80,70 

MADD 760 

70 BEL=B (IJB) 

MADD 7*?0 

80 R (IJR)=AEL+BEL 

MADD 780 

90 CONTINUE 

MAT® 790 

RETURN 

MADD 800 


MADD 810 

AID MATRICES FOR OTHER CASES 

MADD 820 


MADD 830 

100 DO 110 1=1, NM 

MADD 840 

110 R(I)=A (I) +B (I) 

MADD 850 

RETURN 

MADD 860 

END 

MAC® 870 
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c 

c 


MPRD 10 
. . .MPRD 20 

c 


MPRD 30 

c 

SUBROUTINE MPRD 

MPRD 40 

c 


MPRD 50 

c 

PURPOSE 

MPRD 60 

c 

MULTIPLY TWO MATRICES TO FORM A RESULTANT MATRIX 

MPRD 70 

c 


MPRD 80 

c 

USAGE 

MPRD 90 

c 

CALL MPRD (A,B, R,N,M,MSA, MSB ,L) 

MPRD 100 

c 


MPRD no 

c 

DESCRIPTION OF PARAMETERS 

MPRD 120 

c 

A - NAME OF FIRST INPUT MATRIX 

MPRD 130 

c 

B - NAME OF SECOND INPUT MATRIX 

MPRD 140 

c 

R - NAME OF OUTPUT MATRIX 

MPRD 150 

c 

N - NUMBER OF ROWS IN A AND R 

MPRD 160 

c 

H - NUMBER OF COLUMNS IN A AND ROWS IN B 

MPRD 170 

c 

MSA - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX A 

MPRD 180 

c 

0 - GENERAL 

MPRD 190 

r 

V-* 

1 - SYMMETRIC 

MPRD 200 

c 

2 - DIAGONAL 

MPRD 210 

c 

MSB - SAME AS MSA EXCEPT FOR MATRIX 3 

MPRD 220 

c 

L - NUMBER OF COLUMNS IN B AND R 

MPRD 230 

c 


MPRD 240 

c 

REMARKS 

MPRD 250 

c 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRICES A OR 

B MPRD 260 

c 

NUMBER OF COLUMNS OF MATRIX A MUST BE EQUAL TO NUMBER OF 

ROWMPRD 2*70 

c 

OF MATRIX B 

MPRD 280 

c 


MPRD 290 

c 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

MPRD 300 

c 

DOC 

MPRD 310 

c 


MPRD 320 

c 

METHOD 

MPRD 330 


C THE M BY L MATRIX B IS PREMULTIPLIED BY THE N BY M MATRIX A MPRD 340 


C 

AND THE RESULT IS STORED 

IN THE N BY L MATRIX R. THIS IS A 

MPRD 350 

C 

ROW INTO COLUMN PRODUCT. 



MPRD 360 

c 

THE FOLLOWING TABLE SHOWS 

THE STORAGE MODE OF THE OUTPUT 

MPRD 3^0 

c 

MATRIX FOR ALL COMBINATIONS OF INPUT .MATRICES 

MPRD 380 

c 

A 

B 

R 

MPRD 390 

c 

GENERAL 

GENERAL 

GENERAL 

MPRD 400 

c 

GENERAL 

SYMMETRIC 

GENERAL 

MPRD 410 

c 

GENERAL 

DIAGONAL 

GENERAL 

MPRD 420 

c 

SYMMETRIC 

GENERAL 

GENERAL 

MPRD 430 

c 

SYMMETRIC 

SYMMETRIC 

GENERAL 

MPRD 440 

c 

SYMMETRIC 

DIAGONAL 

.GENERAL 

MPRD 450 

c 

DIAGONAL 

GENERAL 

GENERAL 

MPRD 460 

c 

DIAGONAL 

SYMMETRIC 

GENERAL 

MPRD 470 

c 

DIAGONAL 

DIAGONAL 

DIAGONAL 

MPRD 480 

c 




MPRD 490 

c 




.MPRD 500 
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SUBROUTINE MPRD ( A , B , R , N ,M , MSA , MSB ,L ) 
DIMENSION A (1) ,B (1) ,R (1) 

SPECIAL CASE FOR DIAGONAL BY DIAGONAL 

MS=MSA*10+MSB 
IF (MS-22) 30,10,30 
10 DO 20 1=1, N 
20 R(I)=A(I)*B(I) 

RETURN 

ALL OTHER CASES 

30 IR=1 

DO 90 K=1 ,L 
DO 90 J=1 ,N 
R(IR)=0 
DO 80 1=1, M 
IF (MS) 40,60,40 
40 CALL LOC(J,I,IA,N,M,MSA) 

CALL LOC(I,K,IB,M,L,MSB) 

IF(IA) 50,80,50 
50 IF (13) 70,80,70 
60 IA=N*(I-1)+J 
IB=M* (K-l)+I 

70 R (IR)=R (IR)+A (IA) *B (IB) 

80 CONTINUE 
90 IR=IR+1 
RETURN 
END 


MPRD 510 
MPRD 520 
MPRD 530 
MPRD 540 
MPRD 550 
MPRD 560 
MPRD 570 
MPRD 580 
MPRD 590 
MPRD 600 
MPRD 610 
MPRD 620 
MPRD 630 
MPRD 640 
MPRD 650 
MPRD 660 
MPRD 67Q 
MPRD 680 
MPRD 690 
MPRD 700 
MPRD 710 
MPRD 720 
MPRD 730 
MPRD 740 
MPRD 750 
MPRD 760 
MPRD 770 
MPRD 780 
MPRD 790 
MPRD 800 
MPRD 310 
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SUBROUTINE TPRD 
PURPOSE 

TRANSPOSE A MATRIX .AND POSTMULTIPLY BY ANOTHER TO FORM 
A RESULTANT MATRIX 

USAGE 

CALL TPRD (A,B,R,N,M,MSA, MSB, L) 

DESCRIPTION OF PARAMETERS 

A - NAME OF FIRST INPUT MATRIX 

B - NAME OF SECOND INPUT MATRIX 

R - NAME OF OUTPUT MATRIX 

N - NUMBER OF ROWS IN A AND B 

M - NUMBER OF COLUMNS IN A AND ROMS IN R 

MSA - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX A 

0 - GENERAL 

1 - SYMMETRIC 

2 - DIAGONAL 

MSB - SAME AS MSA EXCEPT FOR MATRIX B 
L - NUMBER OF COLUMNS IN B AND R 

REMARKS 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRICES A OR B 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
IOC 


METHOD 

MATRIX TRANSPOSE OF A IS NOT ACTUALLY CALCULATED. INSTEAD, 
ELEMENTS IN MATRIX A ARE TAKEN COLUMNWISE RATHER THAN 
ROWWISE FOR MULTIPLICATION BY MATRIX B. 

THE FOLLOWING TABLE SHOWS THE STORAGE MODE OF THE OUTPUT 
MATRIX FOR ALL COMBINATIONS OF INPUT MATRICES 


C 

A 

B 

R 

TPRD 3S0 

c 

GENERAL 

GENERAL 

GENERAL 

TPRD 390 

c 

GENERAL 

SYMMETRIC 

GENERAL 

TPRD 400 

c 

GENERAL 

DIAGONAL 

GENERAL 

TPRD 410 

c 

SYMMETRIC 

GENERAL 

GENERAL 

TPRD 420 

c 

SYMMETRIC 

SYMMETRIC 

GENERAL 

TPRD 430 

c 

SYMMETRIC 

DIAGONAL 

GENERAL 

TPRD 440 

c 

DIAGONAL 

GENERAL 

GENERAL 

TPRD 450 

c 

DIAGONAL 

SYMMETRIC 

GENERAL 

TPRD 460 

c 

DIAGONAL 

DIAGONAL 

DIAGONAL 

TPRD 470 


TPRD 480 
.TPRD 490 
TPRD 500 
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SUBROUTINE TPRD(A,B,R,N,M,MSA,MSB,L) 

TPRD 510 

DIMENSION A(l) ,B(1) ,R(1) 

TPRD 520 


TPRD 530 

SPECIAL CASE FOR DIAGONAL BY DIAGONAL 

TPRD 540 


TPRD 550 

MS=MSA*10+MSB 

TPRD 560 

IF (MS-22) 30,10,30 

TPRD 570 

10 DO 20 1=1, N 

TPRD 580 

20 R(I)=A(I)*B(I) 

TPRD 590 

RETURN 

TPRD 600 


TPRD 610 

MULTIPLY TRANSPOSE OF A BY B 

TPRD 620 

“ 

TPRD 630 

30 IR-1 

TPRD 640 

DO 90 K=1,L 

TPRD 650 

DO 90 J=1,M 

TPRD 660 

R(IR)=0.0 

TPRD 6 7 0 

DO 80 1=1, N 

TPRD 680 

IF (MS) 40,60,40 

TPRD 690 

40 CALL LOC (I ,J,IA,N,M,MSA) 

TPRD 700 

CALL DOC (I ,K,IB,N,L,MSB) 

TPRD 710 

IF(IA) 50,80,50 

TPRD 720 

50 IF (IB) 70,80,70 

TPRD 730 

60 IA=N*(J-1)+I 

TPRD 740 

IB=N* (K-D+I 

TPRD 750 

70 R(IR)=R(IR)+A(IA)*B(IB) 

TPRD 760 

80 CONTINUE 

TPRD 770 

90 IR=IR+1 

TPRD 780 

RETURN 

TPRD 790 

END 

TPRD 800 


origin^ 
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c 



IOC 

10 

c 



,LOC 

20 

c 



IOC 

30 

c 

SUBROUTINE LOC 


IOC 

40 

Q 



IOC 

50 

c 

PURPOSE 


LCC 

60 

c 

COMPUTE A VECTOR SUBSCRIPT FOR AN ELEMENT IN A MATRIX 

OF 

LOC 

70 

c 

SPECIFIED STORAGE MODE 


LOC 

80 

c 



IOC 

90 

c 

USAGE 


LOC 

100 

c 

CALL LOC (I ,J,IR,N,M,MS) 


IOC 

110 

c 



IOC 

120 

c 

DESCRIPTION OF PARAMETERS 


LOC 

130 

c 

I - ROW NUMBER OF ELEMENT 


IOC 

140 

c 

J - COLUMN NUM3ER OF ELEMENT 


LOC 

150 

c 

IR - RESULTANT VECTOR SUBSCRIPT 


LOC 

160 

r 

N - NUMBER OF ROWS IN MATRIX 


LOC 

170 

c 

M - NUMBER OF COLUMNS IN MATRIX 


LOC 

180 

c 

MS - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX 


LOC 

190 

r 

0 - GENERAL 


LOC 

200 

c 

1 - SYMMETRIC 


IOC 

210 

C 

2 - DIAGONAL 


LOC 

220 

c 



IOC 

230 

c 

REMARKS 


LOC 

240 

c 

NONE 


IOC 

250 

c 



IOC 

260 

c 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 


LOC 

270 

c 

NONE 


LOC 

280 

c 



LOC 

290 

c 

METHOD 


LOC 

300 

c 

MS=0 SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*M ELEMENTS 

LOC 

310 

c 

IN STORAGE (GENERAL MATRIX) 


LOC 

320 

n 

V. 

MS=1 SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*(N+l)/2 IN 

IOC 

330 

c 

STORAGE (UPPER TRIANGLE OF SYMMETRIC MATRIX) . 

IF 

IOC 

340 

c 

ELEMENT IS IN LCMER TRIANGULAR PORTION, SUBSCRIPT IS 

LOC 

350 

c 

CORRESPONDING ELEMENT IN UPPER TRIANGLE. 


LOC 

360- 

c 

MS=2 SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N ELEMENTS 

LOC 

370 

c 

IN STORAGE (DIAGONAL ELEMENTS OF DIAGONAL MATRIX) . 

LOC 

380 

c 

IF ELEMENT IS NOT ON DIAGONAL (AND THEREFORE NOT IN 

IOC 

390 

c 

STORAGE) , IR IS SET TO ZERO. 


LOC 

400 

c 



LOC 

410 

c 



.IOC 

420 

c 

■ 


LOC 

430 


SUBROUTINE LOG (I ,J,IR,N,M,MS) 


LOC 

440 

c 



IOC 

450 


IX=I 


LOC 

460 


JX=J 


LOC 

470 


IF (MS-1) 10,20,30 


LOC 

480 


10 1RX^1*{JX-1)+IX 


DOC 

490 


GO TO 36 


LOC 

500 


20 IF(IX-JX) 22,24,24 


LOC 

510 


22 IRX=IX+(JX*JX-JX)/2 


IOC 

520 


GO TO 36 


IOC 

530 


24 IRX=JX+(IX*IX-IX)/2 


IOC 

540 


GO TO 36 


LOC 

550 


30 IRX=0 


LOC 

560 


IF(IX-JX) 36,32,36 


LOC 

570 


32 IRX=IX 


LOC 

580 


36 IR=IRX 


IOC 

590 


RETURN 


LOC 

600 


END 


LOC 

610 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

w 

c 

r 

n 


MSUB 10 
MSUB 20 
MSUB 30 


SUBROUTINE MSUB 
PURPOSE 

SUBTRACT TWO MATRICES ELEMENT BY ELEMENT TO FORM RESULTANT 
. MATRIX 

USAGE 

CALL MSU3<A,B,R,N,M,MSA,MSB) 

DESCRIPTION OF PARAMETERS 
A - NAME OF INPUT MATRIX 
B - NAME OF INPUT MATRIX 
R - NAME OF OUTFUT MATRIX 
N - NUMBER OF RONS IN A,B,R 
M - NUMBER OF COLUMNS IN A,B,R 

MSA - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX A 

0 - GENERAL 

1 - SYMMETRIC 

2 - DIAGONAL 

MSB " SAME AS MSA EXCEPT FOR MATRIX B 

REMARKS 

NONE 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
LOG 


METHOD 

STRUCTURE OF OUTPUT MATRIX IS FIRST DETERMINED. SUBTRACTION 
OF MATRIX B ELEMENTS FROM CORRESPONDING MATRIX A ELEMENTS 
IS THEN PERFORMED. 

THE FOLLOWING TABLE SHOWS THE STORAGE MODE OF THE OUTPUT 
MATRIX FOR ALL COMBINATIONS OF INPUT MATRICES 

A B R 


MSUB 360 

Mcrra nn 


GENERAL GENERAL GENERAL MSUB 380 

GENERAL SYMMETRIC GENERAL MSUB 390 

GENERAL DIAGONAL GENERAL MSUB 400 

SYMMETRIC GENERAL GENERAL MSUB 410 

SYMMETRIC SYMMETRIC SYMMETRIC MSUB 420 

SYMMETRIC DIAGONAL SYMMETRIC MSUB 430 

DIAGONAL GENERAL GENERAL MSUB 440 

DIAGONAL SYMMETRIC SYMMETRIC MSUB 450 

DIAGONAL DIAGONAL DIAGONAL MSUB 460 

NSU3 4~?0 

' MSUB 480 

MSUB 490 

SUBROUTINE MSUB (A,3,R r N,M,MSA,MSB) MSUB 500 
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DIMENSION A (1) ,B(I) f R(l) 

DETERMINE STORAGE MODE OF OUTPUT MATRIX 

IF (M3A-MSB) 7,5,7 
5 CALL LOC(N,M,NM‘,N,M,MSA) 

GO TO 100 
7 MTEST=MSA*MSB 
MSR=0 

IF(MTEST) 20,20,10 
10 MSR=1 

20 IF(MTEST~2) 35,35,30 
30 MSR=2 

LOCATE ELEMENTS AND PERFORM SUBTRACTION 

35 DO 90 J=1,M 
DO 90 1=1, N 

CALL LOC (I ,J, IJR,N,M,MSR) 

IF (IJR) 40,90,40 
40 CALL LOC (I ,J,IJA,N,M,MSA) 

AEL=0.0 

IF (UA) 50,60,50 
50 AEL=A(IJA) 

60 CALL LOC (I ,J,IJB,N,M,MSB) 

BEL=0 . 0 

IF(IJB) 70,80,70 
70 BEL=S(IJS) 

80 R(IJR)=AEL-5EL 
90 CONTINUE 
RETURN 

SUBTRACT MATRICES FOR OTHER CASES 

100 DO 110 1=1 ,NM 
110 R (I)=A(I)-B (I) 

RETURN 

END 


MSUB 510 
MSUB 520 
MSUB 530 
MSUB 540 
MSUB 550 
MSUB 560 
MSUB 570 
MSUB 580 
MSUB 590 
MSUB 600 
MSUB 610 
MSUB 620 
MSUB 630 
MSUB 640 
MSUB 650 
MSUB 660 
MSUB 670 
MSU3 680 
MSUB 690 
MSUB 700 
MSUB 710 
MSUB 720 
MSUB 730 
MSUB 740 
MSUB 750 
MSUB 760 
MSUB 770 
MSUB 780 
MSUB 790 
MSUB 800 
MSUB 810 
MSUB 820 
MSUB 830 
MSUB 840 
MSUB 850 
MSUB 860 
MSUB 870 
MSUB 880 



EXPERIMENTAL DESIGN FOR OPTIMAL PASS SELECTIONS 


The objective of this task is to develop and demonstrate feature 
selection programs for the purpose of selecting a priori statistically 
optimum subsets of LANDSAT acquisitions for analysis to separate wheat 
from nonwheat when given an adequate sample of labelled wheat and nonwheat 
LACIE seqment pixel data. Throughout the remainder of this report we will 
identi'fy a LANDSAT acquisition, not by its Julian date, but by a sui table 
(e.g. Robinson) wheat growth stage index. Denote the distinct values of 
this index corresponding to growth stages, by 



for winter wheat and spring wheat respectively. T he problem is as follows. 
Select the subsets of size 1, 2 , 3 and 4 of X and X c which maximize the 
separations of spectral measurement of wheat and nonwheat. In the spring 
wheat areas we wish to separate small grains from other crops. 

The Data Set 

The data set will consist of a number of blind sites, preferably 
from more than one growing season. Multiple acquisitions will be required, 
however, we will not require every possible acquisition from every segment.' 
Each acquisition will be assigned the appropriate wheat growth stage index. 

If a segment has multiple acquisitions with the same index, only one (selected 
at random) will be retained. 



The data can now be arranged as follows: 



The V"s denote an acquisition. A "blank" will denote no acquisition 
available. . The data will consist of labelled data as statistical signa- 
tures for wheat/nonwheat categories. We recommend that the "Dots" be " 
used in order to avoid dependence on signatures generated from clustering' 
in CAMS. Furthermore, we know that all signatures will not be available 
for all acquisitions. 
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Measure of Separation and Procedure 


The criteria (figure of merit) used will depend on they type of data 
available. 'The two cases arc: 

1 1 Dots : Set s^{j) -jw-j, •••> w a » °-j » •••» c> b ^ be the set of 

labelled dots for segment^ at crop stage tj. Let Dg(j) be a 
discriminant rule based on the prototypes s ^(tj). The functional 
form of Djjj { j ) wi 1 1 be fixed over all segments ^ and crop stages 
however, the coefficients are estimated for each segment and crop 
stage separately. 

Let CjjU) be the proportion of prototypes misclassi fied by D^(j). 
If n^(j) is the number of prototypes of s^{j), then the figure of 
merit for crop stage t. is defined to be: 

J 


/u) « 


n / J ") 


where A(j) is the index set of all segments having acquisitions 

for crop stage t j , and N(j) -^n g (j ) . That is to say, ^(j) is 

At A<v> 

the misclassification rate over all segments for stage t.. We 

J 

define the best c rop stage for discriminating between wheat and 
nonwheat (small grains and other in spring wheat areas) to be that 
crop stage t. for which ^(j) is a minimum. In order to choose the 

J 

best pair of crop stages, define s^{t.j , t^) to be the set of two pass 
prototypes (8 X 1 vectors) and D^(i,j) the chosen discriminant rules. 
Now A(i,j) is the index set of segments having acquisition at crop 
stages ti. and t. (note that A(i,j) A( i ) . Now 


/(i >J ) 


N jt (i,j)C^(i,j) 


N(i ,j) 





where the definitions of n^(i,j), C^(i,j), and N(i,j) are extended in the 
obvious way. The best pair of crop stages is determined by the maximum 
value of j>. In general, we minimize 



.... \) = 



» • ■ • » ^ ^ ^ 1 » * * * © 


linn >h) n(1 i V 


to choose the best k crop stages. We are obviously constrained by the 
availability of a sufficient number of segments with acquisition at crop 
stages t.-j , t^ for all such subsets of size k of I w or I $ . The 
size of A(i-j, i^) must be greater than two for all i-j , 1^. 

(Interesting values of k are 1, 2, 3, 4.) This constraint will most likely 
require use of aVI_ blind sites in the s.tudy. 

The next is„sue is the definition of the discriminant rules D^(1.|, •«*> ij 
We propose two such rules. Since we have only dots, a -nonparametric discrimi 
nant rule is recommended (although statistical attributes of the categories 
or subcategories might be constructed^ . 


Rule 1 . For each segment find the hyperplane (which may or 

may not separate the two categories of prototypes) determined by 

some non-pa rametric linear discriminate function (LDF) technique. 

Rule 2 . Use the same LDF algorithm to determine a quadratic 

discriminant function. Specifically, transform each acquisition to 

TACAP brightness (b) - greenness (G) space, , “ ' , , G^. . 

2 2 

Construct duimiy variables , G. , for each acquisition so 

that (together with B- and G. ) five variables are defined per 

1 l « 

acqusition for application of the LDF procedure. This will yield 
a quadratic discriminant function based on the B-G coordinate system. 



DATA REQUIREMENTS FOR~ LIMITED TEST 'OF EXPERIMENTAL 
DESIGN FOR OPTIMAL PASS SELECTION . 

The University of Houston will perform a 
limited test of the optimal pass selection experimental 
design, described earlier in this report, in 1979 
with the following data requirements. 


At least 10 blind sites in winter wheat 
strata, which have 4-7 acquisitions, will 
be designated. Labeled dots along with 
NASA specified crop indices for each segment 
and acquisition will be supplied to UH in 
9 track 800 BPI EDCDIC tapes. 

We will perform the limited test and document test 
procedures within the scope of the 1979 contract. 



1. Introduction 


Feature Selection for Best Mean Square 
Approximation of Class Densities . 


The purpose of this note is to describe a general mean square approach 

to linear feature selection which connects certain generalized Fisher criteria 

in discriminant analysis with a measure of pattern class separation intro- 
( 3 ) 

duced by Devijver 7 . The former are typical of those criteria which 
utilize only low order Information about the pattern class distributions, 
while the latter requires' that the class distributions be known, or at 
least accurately estimated. 

Let X denote a random vector in real n-space R n which arises from one 

of m pattern classes , IT having known prior probabilities 

m 

where > 0 and E a. a 1, Let F.(x) denote the j — class 
i=l J m 

conditional distribution function of X and let F(x)> = E a, F (x) denote 

i-1 , 

the mixture distribution. For a given measureable transformation T:R n R 
let G^ (y, T) and G(y, T) denote, respectively, the j — class conditional 
distribution and mixture distribution of the random variable Y = TX. We 
let f (x) (resp. g^(y, T) ) denote the class conditional densities of X 
(resp. Y) with respect to their corresponding mixture distributions; i.e. , 


dF. 

f j “ dF 1 aTld S j ( *’ T) 


dG.(.’ T). 
dG(. , T) 


We will restrict our attention to the set of linear transformations T of 

rank k, and assume that each pattern class Ik has a mean vu and positive de- 
ni 

^finite covariance matrix 9 U . Let p = E a,p . and let 

i . , i i 

i=l 



2 . 



in 



Vi 


- u) (t^ - y) T 


m 

s « m .\ °i a i 

i=l 


and 


S = 


S w + S B 


denote the between class scatter matrix, the average within class scatter 
matrix, and the total scatter matrix respectively. 

A number of interesting feature selection criteria can be formulated 
using only the parameters y, y^, (2 , S, S^, S B ; e.g., the criteria proposed 
by Kittler and Young Foley and Samraon^\ Fukanaga and Koontz^\ 

and the discrete analogue of the modified Karhunen-Loeve expansion of Chlen 
and Fu^-. The modified K-L expansion minimizes an entropy function, and 
also best represents the pattern vector X in an overall least squares 
sense; however, its value for discrimination has been questioned by several 
authors (see Kittler^). Fukanaga^ considers several criteria of the 
generalized Fisher type, including 


J k (T). « tr(T T S w T) _1 (T T S B T). 


Thus, according to this criterion, the best k x n matrix T of rank k is one 
which maximizes J^(T). The solution is any T which is row equivalent to a 
k x n matrix whose rows are linearly independent principal eigenvectors 
(i.e., corresponding to the largest eigenvalues) of S''" Sg. We also consider 


a modification 



3 . 


J k (T) = cr (T T ST) 1 (T T S b T) 
which admits the same maximizing T. 

The Bayesian distance corresponding to the pattern classes 11^ 

( 3 ) 

as defined by Devijver , is 


m 2 " J 2 

B ” £ af E [f . (X) ] 

n ±=i 1 1 

= ” a 2 I f,(x) 2 dF(x). 

i-X 1 R n / 

Its transformed value is 

m 9 o' 

B (T) = l a J g.(y, T) Z dG(y, T) . 

1=1 R 

Devijver proves a number of interesting inequalities relating B^ to the 

Bayes probability of misclassif ication, the Bhattacharrya coefficient, and 

other measures of class separation. In addition, he notes that Cover and 
<2} 

Hart have shown that 1-B is the asymptotic error rate of the nearest 

n 

neighbor classifier. 

2. Mean Square Optimality of Bayesian Distance 

For a given k x n matrix T of rank k, let L^(T) denote the set of all 

k 1 r 2 

measureable functions : R R such that j^k cp(y) dG(y, T) < 00 and let 
0^ be a given closed linear subspace of t^T). Our general approach to 
linear feature selection is to choose that T, y if possible, which minimizes 



4 . 


m ,, 

R(T) * E B min / [cp (Tx) - f (x) ~T dF(x)-, 

i“l * , 

where the 3^ are positive weights* That is, we attempt to find a T which 

produces a set of approximations 9 (Tx) to the class densities f^(x) which 

is best in an overall mean square sense. Given such approximations we may 

classify observations of X according to the pseudo-Bayes rule: decide that 

X is from class IL if a^<p^(TX) > a_. cp ^ (TX) for each j ^ i. Since we are 

interested in classification accuracy, it seems appropriate to choose weights 

3^ which reflect the relative importance of the classes in the mixture 

? 

distribution; e.g., 3, = a, for all i or 8 = a, for all i. For the 

x i i i 

2 

remainder of this section we choose 3^ = andX^, = I^CT). 

2 

Proposition 1 : For 3 i - a^, i =* 1, * . . , m and « L^CT) for each T, 

R(T) = - B fc (T). 

Proof : Observe that g^(y, T) e L^(T), since it is bounded by a^. Moreover, 

for each <pe L^(T), 


/ p(Tx)[g.(Tx, T) - f . (x) 3 dF(x) 

R n 3 1 1 

= / cp(Tx) g (Tx, T) dF(x) - / <P(Tx) dF,(x) 

R n 1 R 

= / . <P( y) g (y, T) dG(y, T) - / . *(y) d G (y, T) 
R fc 1 R 


- 0 . 



Therefore, 


min / t9(Tx) - f, (x)] 2 dF(x) 
9eL 2 (T) R n 


= / n Cg i (Tx, T) - f ± (x)] 2 dF (x) 

= f f (x) 2 dF (x) - / g.Cy, T) 2 dG(y, T) . 
R n 1 R 

2 

The assertion of the proposition follows on multiplying by cu and summing 
over i. 

We may summarize by saying that if there exists a k * n matrix Tq of 
rank k which maximizes then the functions g^d^x, T^) , g^CT^x, T^) 

constitute the best mean square approximation to the class densities 
f^(x) , . f ffi (x) attainable through a linear compression of the data into 
k dimensions. Since B^(QT) = B^d) for each nonsingular k * k matrix Q and 
each k * n matrix T of rank k, B^d) has a maximum if and only if it has a 
maximum on the compact set {T j TT = 1^.^}. particular, if B^d) is 
continuous, it has a maximum. 

3. Best Linear Approximation of Class Densities 

T 

In this section we let C^, be the set of functions 9(y) = w + b y, where 

w is a real number and b c R . For simplicity, we use the notation 

<p(y) = a^v(y), where a = (^) (. and v(y) = (~) c R^ + ^. For given T, 

w ± 

a. ** Cr — ) minimizes 
i 
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/ Ca T v(Tx) - f . (x) ] dF(x) 

V 1 


a [/ v(Tx)v(Tx) dF (x) ]a 
R 


-2a J v(Tx) dF (x) + / f,(xV 

V 1 R n 1 


= a 


T/ 1 


, Hi 


U T T T \ _ n /I I ,.T m T v 

I twt t / 


] a -2(1 | pjT)a 


if and only if 


+ / f., (x) dF(x) 

R n x 


/i|; t t Y " 

V™| TOT T A b 


_ 1 _ 

TV* 


where W = E [XX^] = S + Solving this system gives 


T t T -1 

w ± = i - y t (tst ) - y) 


and 


b, = 5 (tst t )'' 1 t(u. : - yy. 

i i 


The corresponding squared error of approximation is 


_ 1 _ (u - y) T T T (Tsr T ) _1 T(y. - v) + / f.(x) z dF(x). 
i i R n i 


dF (x) 



7 . 


Therefore, the criterion to be minimized is 

m T 

R(T) = - Z & (pi - m) T T T (TST T ) -1 T(M. - u) 
i=l 1 1 


+ terms independent of T. 
That is, we want to maximize 


R(T) = trace (TST^) ^ XS T*' , 

B 

where 

a tn 

S B = 2 3.(M. - M)(M. - y) T . 

i-1 1 x 1 

The solution is X = QTq, where is a k * n matrix whose rows are linearly 
independent principal eigenvalues of , S and Q is an arbitrary nonsingu- 
lar k x k matrix. In particular , for 0^ — we obtain the same solution, 
given by Fukanaga's criterion, 

trace (TS T T ) _1 (TS^T 1 ) . 

W B 


4. Concluding Remarks 

An equivalent set of criteria for feature selection are expressions 
such as 


_ ra _ 

R(T) = Z gi min / [9(Tx) - d. f.(x)J 2 dF(x) 

i=l feC T R n 11 
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in which the posterior probabilities a^f^(x) of the classes are approximated. 

— — 2 
If each Q is chosen to be 1 R(T) is the same as R(T) with 3^ « , 

This, together with Proposition 1 and the relationship between Bayesian 

distance and the probability of error, seems to indicate that the choice 
2 

3 *= is a good one. 

In some cases it may be numerically feasible to use B^(T) as a feature 
selection criterion when assumptions about the parametric form of the class 
distributions are made. For example, if each class distribution F^ is 
multivariate normal, then B^(T) reduces to an expression which is continu- 
ously differentiable in T and which, moreover, can be approximated by 
sample averages over an unlabeled sample from the mixture distribution* 

Thus descent algorithms might be successfully employed in maximizing B^(T). 

Finally, we remark that there is no reason in principle why B (T) 
cannot be regarded as a criterion for nonlinear feature extraction* Indeed, 
Proposition 1 remains true when T is any measureable transformation from 
R n onto R k . 
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A CONVERGENCE CRITERICN 


FOR OPTIMAL LINEAR COMBINATION 


PROCEDURES 


I. Introduction 

In nil that follow:; 9^ will denote the set of n::n unitary matrices, 
i/e will v e' : ui? £4 with the Euclidean topolo-y and recall that ^ is a 
compact cet. The subset^ of %i consirtinq of those matrices of tin 
fora li = I -2 xx 1 with lxll~l, x£ R n (t!ic 'lourcholder transformations) 
i3 a compact subset of *t4 . 

If V :V->Hculs(-;r) ir continuous v.e will define ^ !aax =1 'ft b ^(U) 
and examine, for certain r.erucncoo |u.jJ in *£4, the convergence 
of the sentiences .{VtUiJli ,! to v .. 1C „ . Sequences of this smture have 
been constructed and studied by Dcccll rad milcyf J], Dccell and 
•iarini re], Decell and ' ayolnr in connection with ' '.n::imizir.£ certain 
class separability functions in cittern classification problems. 
Throughout this n ocr^, % and J/ will be as described above. 

Definition 1. Y will be c -lied ^-eloped rovided U^Zl t\n&Yfcj* ^ 

imply there exists some H (dependent on U) such that ^(U)^y(HU) -V ov . 

Definition 2. A somcnce I Vt.) .°2-, in *UL ..ill bo called V-convcr~ent 
provided converges. 

Definition 3. A senue *,co in'U.ill be called a V -householder 

sequence irovidea HeJf and i an integer ir.olv 


(1) y / (u i ) ^ V(u irt ) 

(2) VClUi) « Wu i+1 ) 


II. Conversance 


«*** 

<* POOR <5^ 


Proposition 1. Each V-houreholdcr sequence ia '/-convergent 

and l|uV(U i ) = V (U) « l.u.bftl^) for so ;n U^whic’i is the limit of 
a subcc-ucr.ce of 

Proof: Let {u^jf^bc n V-h vuccholder so -nonce. Since, by defini- 
tion, is ’.onotone incrc-oi v: sc a uoo.ee bounded by 

it must converge to l.\i»b. y(U^). Since ^iv compact, some subse- 
quence ^UL j \r~\ converge to U ell. "ore over, the continuity of 

Y insures tJvt lim^U^ _) = V(U)» that is, {V(U^ is a conver- 
gent subsequence of )V(Ui)>^ a^d the convcr — cnee of the latter 



insures the concl ueion of the proposition. This completes the 
proof. 

Proposition 2. Each V -Householder eenuehcc i -converges to V na ,, 
if and orily if V i c J~f -sloped. 

Proof?. Suppose each f -Householder sequence t -converges to 
If Vis not //'-sloped then there is so e U^^for which. V (U ) x 
and V(HU) - V(U) for each II in//. The ee-uence fu i / i ‘T 1 v ' ith TJ i =tJ 
for each i is cle .rly V -Householder sequence for vhich 

lir.iV(U. )’ * a eonfcr Uicti on. 

^ v i y ' max 7 

If Vis (//-sloped '>nc< • U^f is a r -Householder sequence then, 

according to Proposition. 1, there is sooe some subsequence 

lu. .1 fZ for which lfin*f / (U 1 ) = |/(U) and U = linlL . How if V(U V 
1 i j' 3-d l 1 0 x j 

there exists H in // such that V(U) ^ ^ . Since both V and 

max 

the napping are continuous it follows that V(U) = limV(U. ) 

<£ V(TIU) = limV(HU. ) 4 lixV(U. ,) =V(U). The latter inequality 

3 X j j 3 

is absurd and hence li -[/(U. ) = V • This coiinletes the uroof. 

i 1 " . 

The next proposition sheds some light on the nature of the 
convergence of V -Householder sequences for. //-sloped V. This re- 
sult is of particular importance in applying algorithms based upon 
construction of V -H.o" r . cholder sequences such as those described 
inf 3 J, f 6] and [7] . 

Proposition 3. If is a VrHoucch'older sequence and V is 

//-sloped then exactly one of the following hold: 

(1) ir strictly monotonia (and convergent toV na , x ); 

(2) for some intrgor k, l.u.b. V(HUi-) - ^(Ui.-) (in which case 

“‘i'W- 

V/e will first show that the hypotheses and ~(l) imply (2j). 


Proof : 



2 


If [ V(Uj- )J ic no -j strictly loonotanfp \th9n there is sorae integer ‘ 
k for which V(ro k ) * Y(U k+1 > = V(U lc ) for.,kch / !I-in//nnd hence • 
l.u.h. V(HU k ) ^ V (U k ) . IIov: if V(U k ) ■£. then for aone H in// 

V(u k ) * V(m, z ) * ^ rir>x . •This, of o’ourge, contradicts 

l.u.b. V(HU k ) i V(u : .) hone* V(O k ) - < • 

The hypothesis and ~(2) imply, for each k, th?,t there exist: 
an ll in // for which ^(U k ) * V(HU ?C ).- this'll,’ that 

^ U k) * V(fIU k ) - V(U k+1 ). The cony.- $o. bf/{ )J ^ to 

^nax ^ ^ a conseouen.ee of Pro )Osition 2-. 
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lit. Anolicetion 

'.Ye will first orove two ’iropositiono before doing the 
application. 


Proposition 4. L_et T:R n ->R ,n be a linear transformation. Then' 
there exists an orthonornnl basis {v^, . . . ,v. n } of R ra satisfying 
^T(v, ) ,T(v • )} = 0 whenever l-i^j-n. 

Proof: For v^R 1,1 define f(v) = I!t(v)I|~. Since f is continuous 

there exists v* with llv - || =1 and f(v) = f(v a ) whenever Ivlf =1. For 
v in R m the partial of f in the direction v at v* is 2<(T(v a ) ,T(v)}. 
If <v,v'} =0 eaad <T(v'),T(v)) * 0 then W.L.O.G. (t(v" ) ,T(v)) * 0 

p p 

from which we can dedr.ee that there exists a,b with a‘+b~'=l and 

for which fjT(av+bv a )H ^ ||T(v a )ll“. This would contradict the choice 

of v* and it follows that ^T(v),T(v a )) = 0 whenever (v,v a )=0. 

Since the collection Cv a J of all vectors v for which <V,V-> =0 is 

an m-1 dimensional subspr.ee' of S m then by an inductive argument we 

can assuae the existence of _ n orthonormal basis {v-,-, . . . of 

[v'J 1 satisfying ^T(v i ) ,?(v^ )) “0 whenever l-i^j-m-l. By lettihg 

v =v" then fv, 1 is v‘ e desired orthonorr.i o l basis of IT 1 , 

m '■1 7 p .- 5 

This comnletes the nroof. 


The collection of h-si tension"! subs o ?.cec of R n will be denoted 
by Let T:R n ->R n be a. linear transformation. For IC injg^ 

define V T (K) = ||T(v^)H <L where [v^,-. ..,v^} is an orthonormal 

basis for K. V,p if well defined since 

Vm ( K ) = X trace {(Tv. )(Tv, ) i } = XL trace {(Tv*) (vJf T )} 

1 i=l 4 1 1 i-1 

= XI tracej(T^T) (v. vt )| = trace f(T^T )( XI v-v?) land since X^v.sV? 

l 1 i J 1 i=i 1 1 i=l x 

represents the projection P. r of R n onto K. For K in lot K -1 

iV n 

* 

denote the orthogonal complement of K in R‘*. He will say that 
K ir Cl • .r.t. T if {'f(x) , T(y )} =1 whenever x € K and y€K^ If 
jl T ( y ) II ^ ||T(x)ll whenever x € II, y e 11x11= lyll =X then v.c will say that 
II is 02 w.r. t. T. 
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Proposition 5. - Let T :.i n ->2 n be a linear transformation, and let 


K o be in jS'n- 


i) If K 0 is Cl '.\v ;r v;.r.t. T then V. t (K)^V t (K 0 ) for all K, 

ii) If V t (K q ) ^ YjOO •? o • • some K in then there exists H 
in for V irh Vrp( K 0 ) * Mj.(H(K 0 )). 


Proof: For any IC,K m there exists by Proposition 4 an ortho- 

normal basis of K for which (wn) ,P K (v.\j)^ =0 

v'henever 1-i^ j-lc. .7.L.O.G. we may assume that P-^ (w^)=0 for 


1-i-p and P^. (w.)=0 for o+l-j-k. 
*-0 J 


Let v^= 


?*• ('■< ) 


Vi- 


for 1-i-n. jv-,, . . . .v^} If r.n orbhonoraal systen 


in iC Q which e:: tends to on orthonorcnl basis jv^, . . . ,-v,. j of K 0 , 


u c 2 2 

For 1-i-p there exists, real numbers a^, setiefyina a^*-bv=l and 
a vector vi in of norm 1 for which w . =a . v. -Ho. v? . Suppose now 


l li ii 


that K is Cl and 02 v;.r. t. T. Then 

||T(v.a)|| 2 = a 2 |T(v i )l| 2 + b?|?(v?)H 2 * 2a i h' i (T(v i ) ,T(v? )) 

= a 2 |T(v i )P+ b 2 ||T(vf)ll 2 * a 2 ||T(v.)S 2 + b 2 «T(v. )l| 2 = Ht ( V± )|| 2 

for 1-i-p, and ||t(w^ )|| ' - l|?(v^)l|‘‘ for p+l-^-lc. It follows from 
this that V T (IC) £ V^(K 0 ) for any IC in^ whenever K Q is Cl and C2 

w.r. t. T. 

Suppose now that V r (K Q ) * %,(K). Then for some index i we 
must have ||T(v.n )j| ^ ||t(v^)||. Select 31 in satisfying H(v^)=w^ 
and.lUv^^v^ whenever 1-j-fc and J*i. Such an H must exist since 
<v;.,v.> = 0 whenever j*i. (v^, . . . ,v i _ 1 ,v.a , v i+1 , . . . , v 1; } is an ortho- 
normal system spnnnina "(K 0 ) and since ||T(v.a)j| llT(v^)l| then 
^ t (K 0 ) ^ ^ t (II(K 0 )). Phis completes the proof. 

For a finite f.- iily *)f multivariate norv'l densities each of 
•whose covariance wtri:: is b!ic identity .'.atrix, let {v^_, . • ♦ , v„jJ be 
the collection of pairwise differences of :een vectors for distinct 
' pe.irc. For a kxn ra fc matrix 75 let V(3) be the 3-induced inter- 
class divergence C 1 3 . 2 :on 



V(B) = [ II trace ] - 


mlc. 


It can be shown that there exist3 a unitary matrix U in ?/for 
which V(B) = V( [i k f Z] U) where £l 1c lz] is the kxn matrix whose 
i* 5 * 1 row equals e?*. For that U we can determine that 

A 


Y<B) = .Z|P k (0(v i ))i| 2 where P k is the projection of R n onto the 

spah of the vectors {e^, . . . ,ejJ. Let V be the mxn matrix whose 
i^* 1 rov/ is v^ for 1-i-ra. Since 

||* k (U(Vi))|| 2 » (u T (e 1 ),v i } 2 +-*+(u a? (e ]c ),v i ) 2 then 

V(B). = ^||V(^(e i ))ll 2 . If we let K Q be the span of {e 1 ,...,e k } 

then Yus) = fy(U^(K 0 )). There fore maximizing interclass diver- 
gence in this pas& can be done by maximizing the function 
U -»Vy(U(K 6 )) on Proposition 5 implies that this function 

is ^-sloped. 
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